/*	Interactive Gravitational Simulator
	Project: IGS Core
	File: Node.cpp
	Author: Mike Bantegui <mbante2@gmail.com>, Hofstra University 
	Copyright (C) 2012 - 2013 Mike Bantegui

	This file is part of the IGS software package for simulating
	the N-Body problem in real-time. This was developed as part
	of a departmental honors thesis project at Hofstra University
	in the Spring 2012 semester.

	IGS is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	IGS is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include "Body.h"
#include "Node.h"
#include "ViewNode.h"
#include <algorithm>

namespace IGS
{

// Default construct a node
Node::Node() : Graviton(), Branches(0), CriticalRadiusSq(1E30)
{
}


// Constructs a interior node with the specified contents
Node::Node(const std::vector<Node *> &branches, const double thetaSq) : Graviton(), 
	Branches(std::count_if(branches.begin(), branches.end(), [] (Node *b) { return b != NULL; }))
{
#ifndef PoolAllocator
	IsLeaf = false;
#endif

	for (unsigned i = 0, k = 0; i < branches.size(); i++)
		if (branches[i] != NULL)
			Branches[k++] = branches[i];
	CalculateMultipole();
	CalculateSize(thetaSq);
}

// Constructs a leaf node with the specified contents
Node::Node(const std::vector<Body *> &bodies, const double thetaSq) :  Graviton(), 
	Branches(bodies.size())
{
#ifndef PoolAllocator
	IsLeaf = true;
#endif

	for (unsigned i = 0; i < bodies.size(); i++)
		Branches[i] = bodies[i];
	CalculateMultipole();
	CalculateSize(thetaSq);
}

// Deletes this node and all it's children
Node::~Node()
{
#ifndef PoolAllocator
	if (!IsLeaf)
		for (unsigned i = 0; i < Branches.size(); i++)
			delete Branches[i];
#endif
}

// Returns a read only version of the Node
ViewNode *Node::AsReadOnly() const
{
	ViewNode *result = new ViewNode(Position, Radius);

	for (unsigned i = 0; i < Branches.size(); i++)
	{
		const Graviton *branch = Branches[i];
		ViewNode *node = branch->AsReadOnly();

		if (node != NULL)
		{
			node->Depth++;
			result->Branches.push_back(std::shared_ptr<ViewNode>(node));
			result->Depth = std::max(result->Depth, node->Depth);
			result->NodeCount += node->NodeCount;
		}
	}

	return result;
}

// Node does not participate in the collision timescales
void Node::CheckTimescale(const double timescaleSquared)
{
}

// Node does not currently participate in pairwise interactions
void Node::Pairwise(Body &other, const double softening)
{
	throw "The Fast Multipole Method is currently not available";
}

// Compute the force of gravity due to this Node on the Body
void Node::Direct(Body &other, const double softening) const
{
	const Vector &dr = this->Position - other.Position;
	const double &r2 = Dot(dr, dr);

	// Monopole approximation
	if (CriticalRadiusSq < r2)
		IGS::Direct(other, *this, softening, dr, r2);

	// Refinement
	else
		for (unsigned i = 0; i < Branches.size(); i++)
			Branches[i]->Direct(other, softening);
}

// Calculate the multipole potential contribution from this node
double Node::MultipolePotential(const Vector &dr, const double r1,
	const double r2) const
{
	double pot = 0;

#if MultipolePotentialTerms >= 1
	double quad = 0;

	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			double q2 = 3 * dr[i] * dr[j];

			if (i == j)
				q2 -= r2;

			quad += Quadrupole[i][j] * q2;
		}
	}

	const double r5 = r2 * r2 * r1;
	pot -= quad / (3 * 2 * r5);
#endif

#if MultipolePotentialTerms >= 2
	double oct = 0;

	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 3; k++)
			{
				double q3 = 5 * dr[i] * dr[j] * dr[k];

				if (i == j) q3 -= r2 * dr[k];
				if (i == k) q3 -= r2 * dr[j];
				if (j == k) q3 -= r2 * dr[i];

				oct += Octupole[i][j][k] * q3;
			}
		}
	}

	const double r7 = r5 * r2;
	pot += oct / (5 * 6 * r7);
#endif

#if MultipolePotentialTerms >= 3
	double r4 = r2 * r2;

	double hexa = 0;
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 3; k++)
			{
				for (int l = 0; l < 3; l++)
				{
					double q4 = 7 * dr[i] * dr[j] * dr[k] * dr[l];

					if (i == j)
					{
						q4 -= dr[k] * dr[l] * r2;
						if (k == l) q4 += r4 / 5;
					}

					if (i == k)
					{
						q4 -= dr[j] * dr[l] * r2;
						if (j == l) q4 += r4  / 5;
					}

					if (i == l)
					{
						q4 -= dr[j] * dr[k] * r2;
						if (j == k) q4 += r4 / 5;
					}

					if (j == k) q4 -= dr[i] * dr[l] * r2;
					if (j == l) q4 -= dr[i] * dr[k] * r2;
					if (k == l) q4 -= dr[i] * dr[j] * r2;

					hexa += Hexadecapole[i][j][k][l] * q4;
				}
			}
		}
	}

	const double r9 = r5 * r4;
	pot -= hexa / (7 * 24 * r9);
#endif

	return pot;
}

// Pass multipole terms up along the tree
void Node::AddMultipole(Node &branch) const
{
#if MultipolePotentialTerms >= 1
	branch.Quadrupole += this->Quadrupole;
#endif
#if MultipolePotentialTerms >= 2
	branch.Octupole += this->Octupole;
#endif
#if MultipolePotentialTerms >= 3
	branch.Hexadecapole += this->Hexadecapole;
#endif
}

// Calculate the multipole expansion
void Node::CalculateMultipole()
{
	Mass = 0;
	Position.Zero();
	Velocity.Zero();

	// Monopole term
	for (unsigned i = 0; i < Branches.size(); i++)
	{
		Graviton &branch = *Branches[i];

		const double &m = branch.Mass;
		const Vector &r = branch.Position;
		const Vector &v = branch.Velocity;
		const Vector &a = branch.Acceleration;

		Mass += m;
		Position += m * r;
		Velocity += m * v;
		Acceleration += m * a;
	}
	Position /= Mass;
	Velocity /= Mass;
	Acceleration /= Mass;

#if MultipolePotentialTerms >= 1
	Quadrupole.Zero();
#endif
#if MultipolePotentialTerms >= 2
	Octupole.Zero();
#endif
#if MultipolePotentialTerms >= 3
	Hexadecapole.Zero();
#endif

	// Compute the multipole terms
	for (unsigned i = 0; i < Branches.size(); i++)
	{
		const Graviton &branch = *Branches[i];

#if MultipolePotentialTerms >= 1
		const Vector &dr = branch.Position - Position;
		const double &r2 = LengthSq(dr);
#endif

#if MultipolePotentialTerms >= 1
		// Quadrupole term
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 3; k++)
			{
				double q2 = 3 * dr[j] * dr[k];

				if (j == k) q2 -= r2;

				Quadrupole[j][k] = branch.Mass * q2;
			}
		}
#endif

#if MultipolePotentialTerms >= 2
		// Octupole term
		for (int j = 0; j < 3; j++)
		{
				for (int k = 0; k < 3; k++)
				{
					for (int l = 0; l < 3; l++)
					{
						double q3 = 5 * dr[j] * dr[k] * dr[l];

						if (j == k) q3 -= r2 * dr[l];
						if (j == l) q3 -= r2 * dr[k];
						if (k == l) q3 -= r2 * dr[j];

						Octupole[j][k][l] = branch.Mass * q3;
					}
				}
		}
#endif

#if MultipolePotentialTerms >= 3
		const double r4 = r2 * r2;

		// Hexadecapole term
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 3; k++)
			{
				for (int l = 0; l < 3; l++)
				{
					for (int m = 0; m < 3; m++)
					{
						double q4 = 7 * dr[j] * dr[k] * dr[l] * dr[m];

						if (j == k)
						{
							q4 -= dr[l] * dr[m] * r2;
							if (l == m) q4 += r4 / 5;
						}

						if (j == l)
						{
							q4 -= dr[k] * dr[m] * r2;
							if (k == m) q4 += r4 / 5;
						}

						if (j == m)
						{
							q4 -= dr[k] * dr[l] * r2;
							if (k == l) q4 += r4 / 5;
						}

						if (k == l) q4 -= dr[j] * dr[m] * r2;
						if (k == m) q4 -= dr[j] * dr[l] * r2;
						if (l == m) q4 -= dr[j] * dr[k] * r2;

						Hexadecapole[j][k][l][m] = branch.Mass * q4;
					}
				}
			}
		}
#endif

		// Pass the terms up the tree
		branch.AddMultipole(*this);
	}

}

// Calculate the size of a node using Ridler's criteria. The critical radius is
// computed as being the maximum of: 1) The cell radius, 2) The radius of each
// component plus the distance from the center of mass to said component. This
// treats each node (and node components) as being spherical in shape. In testing,
// it was found that this criteria provides a more accurate critical radius than the
// diagonal length of the rectangle bounding box criteria. We still use the same
// acceptance criteria (distance < radius / theta) but use the more accurate Ridler
// estimator for the cell radius
void Node::CalculateSize(const double thetaSq)
{
	if (Branches.size() > 1)
	{
		for (size_t i = 0; i < Branches.size(); i++)
		{
			Graviton &branch = *Branches[i];

			const Vector &dr = branch.Position - this->Position;
			const double &r2 = LengthSq(dr);

			// Find maximal radius sphere that encloses all the branches
			Radius = std::max(sqrt(r2) + branch.Radius, Radius);
		}

		CriticalRadiusSq = Radius * Radius / thetaSq;
	}

	// One body present -- force direct evaluation
	else
		CriticalRadiusSq = 1E30;
}

}